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ABSTRACT 

We make a detailed comparison of a recent simulation of the formation of a Milky Way- 
like galaxy in a cosmological context with recent observations of how the structural and 
kinematic properties of the Galactic disk vary with stellar abundance. The simulations 
are part of the Making Galaxies in a Cosmological Context (MaGICC) project, a series 
of fully cosmological smoothed particle hydrodynamics (SPH) simulations that use 
early stellar and supernova feedback to match the stellar mass-halo mass relationship 
found in abundance matching studies. The simulated galaxy has a total stellar mass 
of 2 x 10 10 M , and a total halo mass of 7 x 10 11 M , resembling the Milky Way in its 
gross properties: it is disk-dominated with a flat rotation curve, and the scale length 
and height of the disk are similar to the Milky Way's. The simulated stellar disk is, 
however, kinematically ~50% hotter than the Milky Way. 

Examining stellar sub-populations within a narrow [Fe/H] and [a/Fe] abundance 
range, "mono-abundance populations" , we find remarkable qualitative agreement be- 
tween the simulations and observations in several respects: a) The a-enhanced, chem- 
ically old stars lie in a rotating thickened distribution with a short (2 kpc) radial scale 
length while the solar a-abundance (chemically young) stars form a thin (500 pc) dis- 
tribution with a wide range of scale lengths that depend on their [Fe/H] metallicity. 
b) At a given [a/Fe], there is an outward metallicity gradient, c) Mono-abundance 
populations exist with a continuous distribution of scale heights (from thin to thick) 
and of contributions to the stellar surface disk mass. The simulation differs in some 
details for the observed Galactic disk structure: a) The a-enrichment varies from solar 
to [0/Fe]~ 0.25 rather than a peak [a/Fe]~ 0.5 in the Milky Way thick disk, b) In the 
simulations, there is a distinct very thick, rotating disk-like component with a scale 
height of 1.4 kpc, a remnant of the early, turbulent formation phase of the simulated 
disk, which is not observed in the Milky Way. 

This broad agreement between simulations and observations allows us to test the 
validity of observational proxies used in the literature: we find in the MaGICC simu- 
lations that mono-abundance populations are good proxies for single age populations 
(< 1 Gyr) for most abundances. 

Key words: galaxies: formation - galaxies: structure - hydrodynamics - methods: 
N-body simulation 



1 INTRODUCTION 

Ever since the discovery that the luminosity distribution of 
* Email: stinson 'at' mpia.de edge-on SO galaxies could be well fit by including a thick and 
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thin disk component (|Bursteinlll979h . evidence has accumu- 
late d for the ubiquity of thick di sk components in galax- 
ies. Ivan der Kruit &; Sear id (|l98lh found that the disk of 
the late type galaxy NGC 891 required three separate com- 
ponents to fit its vertical profile, a bulge, thick and thin 
disk. Number counts of dwarf stars in the Solar Neighbor- 
hood found that the Milky Way disk also was w ell fit us- 
ing two distinct scale heights (jGilmore fe R eid 1983). Subse- 
quent observations found that every observed disk galaxy are 
characterized by more than one sing le vertical scale height, 
but are well fit by two components (Ivan der Kruit &; Sear id 
19821: iDalcanton fe Bernstein! l200l 120021: ISeth et al] 12005: 



Yoachim fc Dalcantonll200d : Ide Jong ^a L 2007). Recently, 
large surveys have allowed more extensive examinations of 
the Mil ky Way disk struc ture. With the Sloan Digital Sky 
Survey, [juric et all (| 20081 ) confirmed that the spatial struc- 
ture of the Galactic disk can be well described with two 
distinct scale heights. 

In addition to the density structure of the stellar disks, 
observations have also examined the chemical abundances 
of the stars that make up different sub-populations in 
the Milky Way. A number of spectroscopic studies have 
found distinct stellar abundances in star s that are kine- 
matically associated with the thick disk. ([M aiewski 1993 
Gilmore et aDll995l: IFuhrmannl [l998l: IChiba fe Beerslb oOO 
Prochaska et al.ll2000l: iBensbv et al.ll20Q3l. l2Q05|: IWvse et al 



iReddv et al.l 1 20061 : IFuhrmannl l20oI 
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Stars in the thick disk components are more a- 
enhanced than kinematically cooler thin disk stars. It is rea- 
sonable to assume that a-enhanced stars are older since they 
formed from gas not yet polluted with the iron rich eject a 
of longer-lived Type la supernovae. 

Spectroscopic surveys have begun to assemble more 
extensive samples and have thus allowed more detailed 
comparisons between stellar abundance, disk structure, 
and kinematics. The Geneva-Copenhagen Survey found the 
largest radial [Fe/H] gradient in stars with solar [a/Fe] 
abundance ([Nordstrom et al.ll2004r ). The Sloan Digital Sky 
Survey found that there is a smooth transition from th e 
metallicity of thick to thin disk stars (|lvezic et al.l l2QQ8h . 
iNavarro et al.l (|201lh used a compilation of stellar spectra 
from the literature to find a separation between kinemati- 
cally distinct thick and thin disk populations when selected 
based on their position in the [a/Fe]- [Fe/H] plane. 

One of the largest survey to date, SEGUE, compiled 
y 50000 spectra from stars in the solar neighborhood 
(|Yannv et al. II 20091 ). IBovv et al.1 (2012a, hereafter B12) stud- 
ied a large sample of these and divided the stars into mono- 
abundance populations (hereafter, MAPs) each of which 
span a small range in both [a/Fe] and [Fe/H]. B12 stud- 
ied each MAP individually to determine how the structural 
and kinematic disk properties evolved as a function of ele- 
ment abundance. B12 made maximum likelihood fits of the 
disk scale heights and lengths for each MAP. B12 found 
that each MAP's structure can be described using a single 
scale length that varies systematically with both [a/Fe] and 
[Fe/H]. The a-enhanced MAPs had taller scale heights and 
shorter scale lengths than MAPs with solar abundance pat- 
terns. This difference implies that the oldest stars formed 
with shorter scale lengths when the disk was young. Addi- 
tionally, the spatial structure of the Galactic disk can be well 
described by choosing MAPs of approximately solar [a/Fe] 



and showing how their scale lengths depend on [Fe/H]: more 
metal-rich MAPs have shorter scale lengths than [Fe/H]- 
poor ones, which reiterates the presence of a significant ra- 
dial [Fe/H] gradient in the thin, young disk stellar popula- 
tion. Finally, the surface density contributions of the scale 
heights from these populations smoothly decreases from thin 
to thick MAPs, which lead to the conclusion that the Milky 
Way is n ot composed of distinct thin and thick disk pop- 
ulations (IBovv et all 2012b). The continuous distribution 
seems to imply that the thick disk populations did not orig- 
inate from one significant calamitous event. 



1.1 Our simulation to compare 

In this paper, we compare a state-of-the-art disk galaxy 
formation simulation that resembles the Milky Way in its 
gross properties to the intricate chemo- dynamical patterns 
observed in the Galactic disk. If our simulations are realis- 
tic, we can explore the validity of the assumptions under- 
lying the interpretation of observations like the assumption 
that mono- abundance populations (MAPs) are comprised of 
roughly "co-eval" stars. 

The simulated galaxy we study here is drawn from the 
MaGICC project, which simulates galaxies constrained to 
match the stellar m ass-halo mass relationship defined in 
iMoster et all (|2012h . The MaGICC simulations use super- 
nova feedback and early (pre-supernova) stellar feedback to 
limit star formation to yield the present-day stellar mass 
prescribed by the abundance matching technique. The con- 
sequence of the early stellar feedback is to delay star forma- 
tion in a Milky Way mass galaxy from the typical burst seen 
in simulations at z = 4 to a much natter star formation his- 
tory that peaks after two gas rich minor mergers that happen 
at z — 1 (Stins on et al ][20l3). This history corresponds to a 
stellar mass to halo mass ratio that evolves similarly to what 
is found using ab undance matching for high redshift lum i- 
nosity functions (jMoster et al.ll2012l : iBehroozi et al.ll2012h . 
Star formation is reduced because early stellar feedback 
maintains gas at 10,000 K and prevents gas from reaching 
the central 2 kpc of the galaxy (| Stinson et al.ll2013h and su- 
pernova feedback ejects the low angular mome ntum material 
that would normally create a massive bulge ([Brook et all 
2012a). The resulting galaxy has an a flat rotation curve 
and an exponential surface brightness profile with a scale 
length of 4 kpc. 

This early stellar feedback has a number of other im- 
pacts on galaxy formation. It drives metal rich outflows 
that create gaseous halos that match observations of OVI 
in the circum-gal actic medium of star forming galaxies 
( Stin son et al.ll2012h . It can also expand the inner dark mat- 
ter density profile of dark matter in galax ies up to nearly L* , 
producing da rk matter density "cores" ([Maccio et alj|2012 ) . 
iBrook et al.l (2012b) compared four simulated dwarf galax- 
ies with observations and found that all simulations resulted 
in good matches to the observed scale lengths, luminosities, 
and gas fractions. 



1.2 Qualitative Structural Evolution 

Simulations can be used to show images of how galaxies 
evolve. Figure [1] shows two representations of the evolution 



© 0000 RAS, MNRAS 000, 000-000 




5 Gyr < Age < 6 Gyr 4 Gyr < Age < 5 Gyr 3 Gyr < Age < 4 Gyr 




A Thick Disk from Stellar Feedback 3 



u 

Q_ 



10 
5 


2 Gyr < Age < 3 Gyr 


1 Gyr < Age < 2 Gyr 


Age < 1 Gyr 



-5 


1 






-10 










-10-5 5 10 


-10-5 5 10 


-10-5 5 10 




x [kpc] 


x [kpc] 


x [kpc] 



Age > 10 Gyr 


8 Gyr < Age < 10 Gy 


6 Gyr < Age < 8 Gyr | 


5 Gyr < Age < 6 Gyr 


4 Gyr < Age < 5 Gyr 


3 Gyr < Age < 4 Gyr 








2 Gyr < Age < 3 Gyr 


1 Gyr < Age < 2 Gyr 


Age < 1 Gyr 






1 


-10-5 5 10 

x [kpc] 


-10-5 5 10 

x [kpc] 


-10-5 5 10 

x [kpc] 



9.0 ^ 



u 

8.5 o- 



8.0 , 



7.5 ~U> 
c 

CD 

u 
ro 

6.5 t 
in 

6.0 jjj 
q3 

5.5 

o 

5.0 ° 



(a) where stars formed 



(b) where stars are at z = (co-eval populations) 



Figure 1. Projections of the stellar density for stars grouped into nine different age bins. The left panel shows how the present-day 
(z = 0) distribution of the various MAPs. The right panel shows the location of the stars when they formed. In both cases, there is a 
trend from round (large scale height, short scale length) to flattened (small scale height, long scale length) distributions, when going 
from old to young sub-populations. 



of the simulation, panel (a) shows where the stars formed, 
while panel (b) shows where those same stars are at z — 0. 
These images give a qualitative sense that the oldest stars 
formed with a short scale length and tall scale height, while 
the youngest stars have a long scale length disk with a 
short scale height. A comparison with the right panels shows 
that the co-eval populations have not evolved significantly 
in their structure other than some slight thickening. Thus, it 
is possible to look at the structural properties of present day 
stellar populations and make inferences about where those 
stars formed. 

In 331 we examine the structural parameters of mono- 
abundance populations (MAPs). Figure [2] shows the z — 
stellar structure for a sample of MAPs. The selection in- 
cludes representatives from the three main MAP families 
we find in our quantitative analysis in old-thick disk, 
intermediate, and young-thin. The upper-left panel repre- 
sents the old-thick population that has a tall scale height, 
short length and enhanced [O/Fe]. The right two panels in 
the upper row represent the intermediate population that 
arises from a significant transition in the evolution of the 
stellar disk. In this population, the stars with the tallest 
scale heights also have the longest scale lengths. The bot- 
tom two rows represent the young-thin population. All of 
these populations have as thin a disk as can be resolved in 
the simulation. The scale lengths, however, vary greatly as 
a function of [Fe/H], creating a significant [Fe/H] gradient 
that is apparent in the youngest, solar [O/Fe] abundance 
populations. 

In what follows, we present a quantitative study of the 
disk structure of the s imulated L* g alaxy, gl536 that was 
originally described in Sti nson et al. (2013), but is briefly 
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Figure 2. Edge-on stellar surface densities for a selection of 
mono-abundance stellar populations (MAPs). 



reviewed in 32 along with the physics used in the simulation. 
In 33J we show how the disk structure of the simulations 
compares with observations. In particular, we compare the 
correlations between scale heights and lengths and chemical 
abundances, [a/Fe] & [Fe/H], and the mass weighted disk 
scale height distribution. 
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2 METHODS 
2.1 Simulations 

We use gl536, a cosmological zoom simulation drawn from 
the McMaster Unbiased Galaxy Simulations (MUGS). The 
simulation starts at z = 99 from a cosmologic ally moti- 
vated mat ter distribution base d on WMAP3 (Sp ergel et al.l 
l2007h : see lStinson et all (|2010h for a complete description of 
the creation of the initial conditions. It includes metallicity- 
dependent gas cooling, star formation, and a detailed chemi- 
cal enrichment model that allow us to make a comprehensive 
comparison with observations of the Milky Way. 

gl536 has a virial mass of 7 x 10 11 M©, a spin parame- 
ter of 0.017, and a last major merger at z = 2.9. The mass 
of the Milky Way dark matter halo is u ncertain, with val- 
ues from 7 x 10 1 * M© to 1.4 x IP* 2 M ([Klypin et alJl2QQ2b 
IXue et al.ir2008l : iBovv et al.ll2012h having been determined 
using a variety of means. gl536 is on the low end of this 
mass range, and thus may not make a perfect comparison 
to the Milky Way. However, a comparison may give an idea 
of which physical processes are important for shaping the 
structure of a galaxy disk. The merger history of gl536 is rel- 
atively quiet, which m ay also be the case for the Milky Way 
(Ha mmer et al.ll2007h . The mergers that do accrete onto the 
galaxy follow prograde orbits. Thus, gl536 maintains a disk 
on a relatively constant spin axis orientation throughout its 
evolution. 

The simulation uses the smoothed particle hydrody- 
namics (SPH) code gasoline (Wads lev et all 2004). Details 
of the physics used in the MaGICC project are detailed 
m btmson et al.1 (|2013h . Briefly, stars are formed from gas 
cooler than T m ax = 1.5 X 10 4 K, and denser than 9.6 cm -3 
according to the K ennicutt Schmidt Law as described in 
IStinson et al. with a star formation efficiency param- 

eter, c*=0.1. The star formation density threshold is then 
set to the maximum density at which gravitational insta- 
bilities can be resolved, / as (n th > 9.3 cm ), where 
M gas — 2.2 x 10 5 Mq and e is the gravitational softening 
(310 pc). 

The star particles are 5 x 10 4 M©, massive enough 
to represent an entire st ellar popu l ation consisting of stars 
with masses given by the Chabriei] (|2003T ) initial mass func- 
tion. 20% of these have masses greater than 8 M© and 
explode as Type II supernovae from 3.5 until 35 Myr af- 
ter the star forms , based on the Padova stellar lifetimes 
([Alongi et al. 1993; Bres san et ai1ll993h . Each supernova in- 
puts Esn = 10 51 ergs of purely thermal energy into the sur- 
rounding gas. This energy would be radiated away before it 
had any dynamica l impact b ecause of the high density of the 
star forming gas (|Kat al 19921 ). Thus, the supernova feedback 
relies on temporarily disabling cooling based on the subgrid 
appro ximation of a blastwave as described in Sti nson et al.l 
(|2006h . 

The stars also chemically enrich the ISM during their 
evolution through the explosions of SNII and SNIa. SNII 
chemical enrichme nt is based on linear fits as a function of 
star mass for the IWooslev fc Weaver! (| 19951 ) model of SNII 
explosions. Linear fits do not suffer from the unrealis tic oxy- 
gen e nrichment that some of the power law fits do ([Gibson! 
120021 ). The oxygen fit used is 

Y Qx = 0.21M* -2M (1) 



and the iron fit is 

Y F e = 0.003M* (2) 

SNIa enrichmen t foll ows the W7 model from 
iThielemann et al.1 (1986), where each SNIa produces 
0.74 M of iron and 0.143 M oxygen. 

The supernovae feedback does not start until 3.5 Myr 
after the first massive star forms. However, nearby molec- 
ular clouds show ev idence of being blown apar t before any 
SNII h ave exploded. IPellegrini et al.l ([201 if ) and lLopez et al.l 
([201 if ) emphasized the energy input from stellar winds and 
UV radiation pressure in 30 Doradus prior to any SNII ex- 
plosions. Thus, in the time period before supernovae start 
exploding, we distribute 10% of the luminosity produced in 
the stellar population (equivalent to the UV luminosity) to 
the surrounding gas without disabling the co oling. Most of 
the en ergy is immediately radiated away, but Sti nson et al.l 
(|2Q13h shows that this early stellar feedback has a signif- 
icant affect on the star formation hist ory of a Milky Way 
mass galaxy and places the halo on the Most er~et al.l ([2012) 
stellar mass-halo mass relationship at z — 0. 

2.2 Stellar sub-populations 

Throughout the paper, we divide star particles into sub- 
populations using either their abundance or their age. We 
refer to these sub-populations as either a MAP (mono- 
abundance population) or a co-eval population, respectively. 
These are defined as follows: 

• Mono- abundance population (MAP): Each population 
is taken from a section of the [0/Fe]-[Fe/H] plane, with 
A[O/Fe]=0.05 and A[Fe/H]=0.1. 

• Co-eval: Stars all formed at the same time. The simula- 
tion is divided into 50 co-eval populations with equal num- 
bers of star particles. Each population contains 11410 stars 
and has a characteristic age-spread of ~ 250 Myr. 

2.3 Fitting the structure of sub-population disks 

For comparison with observations, disk structural parame- 
ters (scale length and height) are fit throughout this paper 
using a maximum likelihood method. Specifically, we fit the 
number density of simulated star particles following B12, 
who fit a number density of G dwarf stars. The radial den- 
sity profile is modeled with an exponential function while 
the vertical density distribution uses a sech 2 (z) function to 
capture the the flattening of the vertical density distribu- 
tion near the disk midplane. The likelihood for finding a 
simulated particle at radius R xy ,i and height \zi\ is 

ln(L) = Ef-^+ln(sech 2 (M)) 
- \n(47rr exp h z (-e rex p (r exp + rmax) + 

- r min (3) 
g rexp (Texp ~\~ Train 

))) 

(tanh(^)-tanh(^)), 

where h z and R exp are the disk structural parameters to 
be fit in the annulus from Tmin = 6 kpc to r macc = 10 kpc 
in radius and z m in — to z m ax = 3 kpc in height. This 
likelihood function is turned into a minimization problem 
by taking its negative. Initial best fit parameters are found 
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using Powell's method for minimization (|Press et al.ll200/t ). 
These parameters are u sed as initial conditions for a n en- 
semble MCMC sampler ([Foreman- Mackev et aj]|2012h . The 
MCMC program samples the posterior distribution. The me- 
dian of the distribution is used as the best fit value because 
it provides better fits than the results from Powell's method. 
The quoted errors are the interval that contains 68% of the 
posterior distribution for the given parameter. 

The profiles for a selection of example MAPs are pre- 
sented in Figure [3] It shows the great variation of profiles 
apparent in the various MAPs, differences that will be de- 
scribed in ^3.11 The models provide a good description of 
the data. 



3 RESULTS 

With the procedure in place for fitting the scale height and 
length of sub-samples of stars in the simulations at z = 0, we 
can now compare the mono-abundance populations (MAPs) 
between the observed Galactic disk and the simulations. We 
further study how the structure of the MAPs compares to 
co-eval populations in 



3.1 Simulated MAP Structure 

Figure [H shows how the disk scale height (left) and length 
(right) vary as a function of [Fe/H] and [O/Fe] for differ- 
ent MAPs in the simulation. There are smooth changes in 
both the scale height and length as a function of abun- 
dance. The tallest scale heights exist in the [Fe/H] -poorest 
and most oxygen-enhanced populations. The scale height 
shrinks as the [Fe/H] -enrichment increases and the oxygen- 
enhancement decreases, reaching a minimum of ~ 200 pc. 

Regarding the scale lengths in the right panel, the 
oxygen-enhanced stellar populations all have comparable 
scale lengths around 2 kpc. In less oxygen-enhanced pop- 
ulations, the scale length generally grows. However, the 
[Fe/H] -enriched stellar populations at solar [O/Fe] have 
scale lengths that are shorter, r exp ~ 1 kpc, than the oxygen- 
enhanced populations. 

The [Fe/H] -enriched populations represent the most 
recent star formation, occurring near the galactic center. 
There is a limited radius inside which the ISM has been suf- 
ficiently enriched to produce such stars. Outside this radius, 
star formation is also ongoing, but with a scale length much 
larger than fits on the colorbar scale. The change in scale 
length is separated by a variation of only 0.2 dex in [Fe/H]. 

The a-enhanced populations show much less varia- 
tion in their scale lengths, staying around 2 kpc. They 
show only moderate variation in their scale heights, which 
gradually decrease with higher [Fe/H] and lower [O/Fe]. 
These are the stars that make up the thick disk of gl536. 
Liu & van de Venl (|2012h found that in the Milky Way, stars 
with a-enhanced abundances share many kinematic proper- 
ties. 

The range in [O/Fe] abundances is not as wide in gl536 
as in the Milky Way (0-0.25 versus 0-0.5). This difference 
may indicate a problem with the exact treatment of chemical 
enrichment in the simulation. 

Figure [5] details the relation between the scale height 
and the scale length for the simulated MAPs (bottom panel) 
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Figure 5. Present-epoch scale height vs scale length for var- 
ious MAPs in the simulations. The stars are split into mono- 
abundance populations, 0.05 wide in [O/Fe] and 0.1 in [Fe/H]. 
The simulated points are coloured by [O/Fe], the observed points 
are all grey. In many cases, the errors are smaller than the sym- 
bol size. The simulated MAPs show some trends. The scale height 
decreases with [O/Fe], while the scale length remains around a 
constant value of 2 kpc. Below a scale height of 600 pc, the scale 
lengths are longer. 



and compares it to the Galactic observations from B12 (top 
panel). The distribution of simulated MAPs shows some 
structure. There is a cluster of a-enhanced points with short 
scale lengths and large scale heights, comparable to the clus- 
ter of a-enhanced points in the top panel of Figure [5] (or 
Figure 5 of B12). 

Figure [5] illustrates an important qualitative difference 
between the simulation and the Galactic disk: the simulated 
disk is ^50% thicker than the Milky Way. The extra thick- 
ness may be due to a difference in mass between gl536 and 
the Milky Way. As mentioned before, gl536 is at the low 
end of the Milky Way total mass estimates and has less stel- 
lar mass than the Milky Way disk. Less mass would make 
the disk potential shallower and thus the disk thicker. Al- 
ternatively, the thicker disk could be a sign that the stellar 
feedback is having too strong an effect on the disk. 

Switching focus to the MAPs with low (solar) [O/Fe], 
Figures 3] & [5] show that they all have scale heights < 500 
pc and a wide range of scale lengths from 1 kpc to large 
values. Figure [3] shows that the low [Fe/H] bins have the 
longest scale lengths, such that the fraction of stars with 
such low [Fe/H] values increases further out in the disk. Fig- 
ure [3] shows that the stellar density profiles are nearly flat 
as a function of radius for these populations. As [Fe/H] in- 
creases, the scale length becomes dramatically shorter. Such 
a scale trend with [Fe/H] is also clear in Figure 5 of B12 and 
points to a metallicity gradient present in recently formed 
stars. 

Our simulation predicts the existence of a compact 
metal-rich stellar population with short scale height and 
length. It is plausible that this population did not show up 
in the analysis of B12 because the SDSS G-dwarf data sets 
do not include stars within 300 pc of the disk midplane or 
sightlines towards the Galactic centre. 

There are also intermediate (0.1-0.15) [O/Fe] stellar 
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Figure 3. The density profiles for a selection of simulated mono-abundance populations. This selection is composed of three columns each 
2 plots across. The left plot of the pair is the vertical fit that uses the sech 2 (z) function while the right column shows the corresponding 
exponential radial scale length fit. The dashed green lines represent the best fit. The vertical error bars on the points are the Poisson 
errors. The top row includes fits from two different [O/Fe] values. The first pair is for [O/Fe] =0.2-0.25, which shows typical profiles of 
the a-enhanced stellar population. The second two pairs have [O/Fe] =0.1-0. 15 and show some variation in scale length, but to a lesser 
degree than the bottom two rows. The plots in the bottom two rows contain a selection of fits with constant [O/Fe] and increasing [Fe/H] 
from left to right. The bottom row is [0/Fe]= 0-0.05, the middle row is [O/Fe] =0.05-0.1. These [O/Fe] values are selected because of the 
large variation in scale length as a function of [Fe/H] as shown in Figure [4] In several cases, the stellar density is nearly flat as a function 
of radius. In these cases, the quoted scale lengths are taken from the bottom 1% of the posterior distribution so represent an extreme 
lower limit of the scale length. 



populations that show a direct correlation between scale 
height and length. Figure |4] shows that this is again related 
to the metallicity gradient of the disk. The most enriched 
MAP ([Fe/H] ~ 0) in this [O/Fe] range has the shortest 
scale length and height of these populations. At progres- 
sively lower [Fe/H] abundances, the scale height and length 
both become longer. This correlation is also apparent in Fig- 
ure 5 from B12, although only in the observed MAPs where 
fits are poor. We will examine this intermediate population 
further in the Discussion 8??. 



3.2 Radial abundance gradients 

To make the [Fe/H] radial gradient clearer and identify the 
best way to search for its signature in observations, Figure 
[6] shows the [Fe/H] gradient in populations divided using 
either age or [O/Fe]. As expected from Figures |4] and [5] 
there is a strong radial metallicity gradient in populations 
that are relatively young, or have solar [O/Fe] values (blue 
symbols), with A[Fe/H]/AR~ 0.04 dex/kpc. In contrast, 
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Figure 6. The metallicity gradient of stellar populations when 
split either by age (left panel) or by ai-enrichment [O/Fe]. Each 
plot shows the mean [Fe/H] in 1.5 kpc radial bins. The radii of 
alternating populations are offset by 750 pc so that the error 
bars can be seen. The vertical error bars represent the spread 
(variance) of the [Fe/H] distribution in each bin. 
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(a) Scale height 
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Figure 4. Dependence of the scale height (left) and length (right) as a function of [Fe/H] and [O/Fe] for different simulated mono- 
abundance stellar populations. Stars are divided into abundance populations with A[O/Fe]=0.05 and A[Fe/H]=0.1. These are the same 
fits as Figure [8] only projected into the abundance plane. 



the old or [O/Fe] -enhanced stars (purple symbols) barely 
have any radial [Fe/H] gradient. Thus, the [O/Fe] -separated 
[Fe/H] gradient is a good proxy for the gradient that appears 
in the youngest populations. 

However, the lack of a significant gradient in older pop- 
ulations makes it difficult to use [Fe/H] gradients to identify 
what causes the disk to thicken. 



3.3 Comparing abundances and ages 

In the simulation, we can examine the age-distribution of the 
stars within each MAP. The left panel of Figure shows the 
mean age for all the MAPs. It shows that the [O/Fe] =0.1- 
0.15 MAPs formed between 5 and 7.5 Gyr ago. In this range 
of [O/Fe] values, the stars that formed the longest ago have 
the longest scale length and height. This is an indication 
that the disk did not strictly form from the inside out. This 
intermediate [O /Fe] population seems to indicate that there 
is a transition epoch as the disk shifts from a short, thick 
distribution to a long, thin distribution. 

The right panel of Figure shows, however, that most 
stars within each MAP did form at similar times, with a typ- 
ical age spread of less than 1 Gyr. This implies that [O/Fe] 
and [Fe/H] jointly are a good proxy for co-eval populations 
across much of the [O/Fe]- [Fe/H] abundance plane. 

The MAP with 0.05<[O/Fe]<0.1 and -0.3<[Fe/H]<- 
0.2 has an exceptionally large dispersion, containing con- 
tributions from stars with a range of 6 Gyr. We overplotted 
the mean [O/Fe] and [Fe/H] values for stars split into age 
bins. These points show that the normal enrichment pattern 
makes an excursion towards higher a-enrichment about 7 
Gyr from the end of the simulation. A look a t the star for- 
mation history of this galaxy in Figure 9 of IStinson et al.l 
(2013) shows that this is the time at which the galaxy un- 
dergoes an increase in star formation. An examination of the 
evolution of the galaxjQ shows that there are two prograde 



http:/ /www. mpia.de/ stinson/magicc/movies/c.ltd.05rp.l/far.mp4 




Figure 8. Present-epoch scale height vs scale length for various 
co-eval populations in the simulations. The stars are split into 50 
equal sized age bins. There is an obvious trend when using the 
age binning. The scale height decreases gradually, while the scale 
length remains around a constant value of 2 kpc. When the scale 
height decreases to 600 pc, the scale length starts to grow. 



minor mergers in rapid succession during this epoch. This 
increase in star formation causes stars to form with higher 
[O/Fe] and then continue on a parallel abundance trajectory. 



3.4 Structure of the Co-eval Populations 

We repeat our analysis of the evolution of structural param- 
eters, but using coeval population rather than MAPs. Figure 
[8] shows the scale height, h z , as a function of scale length, 
r exp for 50 different co-eval populations between Tmin — 6 
kpc and r m ax — 10 kpc. The structural evolution of the 
co-eval populations is smoother and more gradual than for 
the MAPs. The oldest stars are in a thick distribution with 
a scale length of 2 kpc. The younger populations have pro- 
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Figure 7. (a) The mean age of the simulated mono-abundance populations; (b) the age dispersion (central 68% of the age distribution) 
for each mono- abundance population. The dots in the right panel are the mean [O/Fe] and [Fe/H] values for stars binned by their age. 
The age dispersion is less than 1 Gyr for nearly every MAP except [O/Fe]=0. 05-0.1, [Fe/H]=-0.3 to -0.2. This abundance corresponds 
to the time when the abundance evolution moves to a trajectory at higher [O/Fe] due to a sudden increase in star formation. 



gressively shorter scale heights, and longer scale lengths. For 
these young stars, the disk scale length begins to increase 
while the scale height becomes more constant. 

Compared to the MAPs in Figure [5] the wide range of 
scale lengths for the recently formed, solar [O/Fe] MAPs dis- 
appear. Instead, the scale lengths of the youngest population 
are 4-6 kpc, an average of the wide range in the MAPs. Fig- 
ure [6] shows that the [Fe/H] gradient is still quite strong in 
the youngest population, even stronger than for the MAPs 
since the MAPs mix populations of different ages. What is 
apparent from a comparison of the structure of the co-eval 
and mono-abundance populations is that the structure of 
MAPs is a mixture of disk evolution and enrichment. The 
enrichment is not uniform. 

Figure [9] shows the fits of the disk scale height and 
length of stars in an annulus 6 to 10 kpc from the galactic 
center, comparable to the plots in Figure [3] These sample 
fits show the gradual transition from a thick disk with a 
short scale length to a thin disk with a long scale length. 

3.5 The Mass- Weighted Scale Height Distribution 

We can calculate how much different sub-populations con- 
tribute to the local surface density at the solar circle. For 
each co-eval population, we calculate its integrated surface 
density using its scale height and length in the density equa- 
tion : 



p(R, z) — poe Rex p sech ' 



2h z 



integrated over \z\ from to 3 kpc: 



Xr©(/i 2 ) = 2p e 1 



/*3 kpc 

/ ; 

J kpc 



sech A 



2h z 



dz 



(4) 



(5) 



We then plot this surface density as a function of the 
characteristic scale heigh of each co-eval population, which 
roughly maps to metallicity and abundance (Figure 0. Fig- 
ure [10] shows that there is a nearly flat distribution of the 



surface density contributed by populations with large scale 
heights up to 800 pc. There is then an absence of surface den- 
sity contributed by populations with scale heights between 
800 and 1200 pc, but then a distinct component with scale 
heights from 1.2 kpc to 1.5 kpc, which corresponds to the 
old, thick disk formed in the simulation. These scale heights 
are far greater than any so far observed in the Milky Way 
and may be due to the low mass of gl536 or the overheat- 
ing from the stellar feedback. The absence of populations 
with scale heights around 1 kpc illustrates how the simula- 
tion transitions suddenly from it thick disk forming phase 
to thin disk. The data from the Milky Way does not show 
any such sudden transition. 



4 DISCUSSION 

Based on the comparison of our simulations with observa- 
tions, we identified three distinct stellar populations. The 
simulation shows clear indications of the evolution from one 
population to the next as the disk evolves and the signature 
of these transitions are also apparent in the observed MAPs. 



4.1 Old thick disk 

The MAPs with [O/Fe] > 0.1 5 all have constant scale lengths 
of rsj 2 kpc in both the simulated galaxy and the Milky Way. 
These MAPs also have the tallest scale heights, though they 
shrink by a factor of two as the disk evolves. The MAPs are 
all older than 8 Gyr and over the 5 Gyr of their evolution, 
their [Fe/H] increases from -1 to -0.5 in both the simulated 
galaxy and the Milky Way. The shrinking scale height is 
steady as [Fe/H] becomes more enriched in the Milky Way. 
While the shrinking scale height is not so clear as a function 
of [Fe/H] in the simulated galaxy, it is apparent as a function 
of age in Figure [8] 
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Figure 9. Scale height and scale length of the simulated co-eval 
stellar populations. The left column is the vertical fit of a sech 2 {z) 
profile while the right column shows the corresponding radial ex- 
ponential scale length fit. The vertical error bars on the points 
are the Poisson errors. 
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Figure 10. Surface density at the solar radius (8 kpc) as a func- 
tion of scale height for coeval populations. The surface densities 
for the points colored by age were calculated using the scale height 
and length of each co-eval population in the density equation in- 
tegrated over \z\ from kpc to 3 kpc. The points are then placed 
according to their scale height. The solid line represents how the 
surface densities of the simulated galaxy are distributed across 
150 pc bins. The dashed line represents the same distribution of 
surface densities for MAPs taken from lBovv et all (2012b). 



4.2 Intermediate population 

There is a general evolution of scale height as an inverse 
function of scale length. There is one series of MAPs that 
does not follow this trend, but instead scale height directly 
follows scale length. These MAPs all have [O/Fe] =0.1-0. 15 
in the simulations and 0.25 in the observations. However, 
this trend does not show up at all in Figure [8] which displays 
scale height as a function of scale length for co-eval popu- 
lations. It shows a strictly inversely proportional evolution. 
Figure [3 shows that the reason these MAPs show the odd 
trend is because they include the widest range of ages. These 
MAPs are filled with stars ~6 Gyr old that formed when the 
abundance evolved to a more [Fe/H] -enriched track. This 
abundance evolution happens due a suddenly increased star 
formation rate. These populations are much more massive 
in the simulation than in the Milky Way. 



4.3 Young thin disk 

At the lowest [O/Fe] enhancements, the disks in both gl536 
and the Milky Way are their thinnest. For gl536, these pop- 
ulations have [0/Fe]< 0.1 and [Fe/H]> -0.5, while in the 
Milky Way, they have [0/Fe]< 0.2 and [Fe/H]> -0.7. The 
scale heights of these populations range from 500 down to 
200 pc in the Milky Way and from 700 down to 300 pc 
in gl536 with the lower [O/Fe] enhanced populations have 
the thinnest disks. The scale lengths vary greatly for these 
MAPs from 2 to > 5 kpc in the Milky Way and from 1 to 
> 10 kpc in gl536. In both cases, these scale lengths vary 
directly as a function of [Fe/H]. All these MAPs are < 3 
Gyr old in the simulation, so the wide variation in [Fe/H] 
represents the metallicity gradient evident in the young thin 
disk. 
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4.4 Differences between simulation and Milky 
Way 

The simulated galaxy has scale heights that are twice those 
observed in the Milky Way The simulation also shows evi- 
dence for a sudden event that makes the disk much thinner 
over a short period of time. It leaves an absence of stars 
that have a scale height of 1 kpc in the distribution of scale 
heights. The Milky Way plot made from MAPs shows no 
similar absence, so it appears that the MW more smoothly 
evolved from thick to thin disk formation. 



5 CONCLUSIONS 

We make a detailed comparison between the stellar struc- 
ture of a disk formed in a cosmological simulation and recent 
observations of the chemical and kinematic structure of the 
Galactic disk. In particular, we compare the systematic vari- 
ations of disk scale height and lengths between simulations 
and observations. At a qualitative level, we find a striking 
similarity. 

Our disk structural analysis is based on both age and 
chemical enrichment. The oldest stellar populations have 
short scale lengths and tall scale heights. When the the disk 
is decomposed according to metal abundance as in obser- 
vations, the thickest populations that have the short scale 
length are those that are a-enhanced. It is clear in the sim- 
ulations that these are the oldest populations. The young 
stellar populations that have solar a-enrichment have the 
shortest scale heights. Since we have divided our popula- 
tions in [O/Fe] and [Fe/H] simultaneously, we can see that 
these young populations have a wide range of scale lengths 
that is strongly dependent on [Fe/H]. This illustrates how 
much chemical enrichment varies at different locations in 
the disk. In particular, the center of the galaxy is the most 
the metal ([Fe/H]) enriched while the outskirts of the disk 
follow a lower [Fe/H] path to solar [O/Fe]. 

The correspondence between the Milky Way and the 
MaGICC simulation add support to scenarios where thick 
disks form kinematically hot. In future work, we will show 
how the morphology of the MaGICC galaxy evolves, but 
since it is a relatively low mass galaxy, it does not have large 
clumps of star formation at high redshift. Rather, the stellar 
feedback heats the galactic gas early in its evolution so that 
it does not settle into any sort of thin disk until z ~ 1. 
The remaining presence of this thickened stellar component 
is what compares best with the thick disk observed in the 
Milky Way. As the galaxy increases in mass the gas settles 
down and forms a thin disk. 
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